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ABSTRACT 



A comparison of two versions of coupled physical/biological models of the 
oceanic mixed layer and upper pycnocline are presented: the first model (NPZ), based on 
Steel (1977), contains components for nutrients, phytoplankton and zooplankton; a 
second model (NPZD) based on equations from Flierl and Davis (1993) contains 
additional components for detritus and light budgets. These models were used to study 
short duration dynamic mixing events and the resulting biological response. This work 
demonstrates the differences between a basic NPZ model and a NPZD model with the 
inclusion of solar radiation. The goal of this study is to determine the effects of detritus 
and light on phytoplankton and zooplankton production on the short time scale. The 
model solutions show that biological modeling is very dependent on physical forcing 
parameters such as turbulent mixing and light attenuation in the water column. The 
effects of detritus are minimal in the short time scale when compared with other 
parameters. 
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I. INTRODUCTION 



In the past several years oceanographers have made significant advancements in 
understanding the physical processes that influence the organisms that reside in the upper 
ocean. This knowledge comes both from advances in modeling techniques and in 
improved sampling methods. 

This study will continue from the work of Doney et al. (1996) and will focus on 
the advances in coupled physical/biological numerical modeling of the ocean mixed layer 
with an emphasis on the biological aspects. The biological effects following short 
duration transient events that result from dramatic changes in mixed layer forcing, such as 
strong winds will be examined. The results from two coupled physical/biological will be 
contrasted. The first model is a mixed layer model which represents the most basic 
biological interactions, with no component for detritus; the second model solves a more 
advanced set of equations with model components for both detritus and light. 

A brief history and overview of biological modeling is given in Chapter D. In 
Chapter IE the two models used here are described. Chapter IV discusses the results of 
these models, and Chapter V will summarize this work and make some recommendations 
for future work. 
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II. BACKGROUND 



Our understanding of the physical, biological and chemical interactions that occur 
in the world’s oceans has grown over the past several decades. This fact, coupled with an 
exponential growth in computer power, has enabled us to make great progress towards 
accurately modeling the physical and biochemical phenomena that are occurring in the 
ocean. This chapter will look at the development of biological modeling starting with 
Steel ini 977 and progressing to current work by Doney (1996) and others. 

Steel (1977) determined the most basic equations for biological interactions in the 
ocean. He pointed out many of the problems in biological modeling, the most obvious 
being the complex and, in some cases, poorly understood interactions among biological 
organisms. The interactions between physical and biological variables is still under 
investigation today and will continue well into the future. Steel’s study provided the 
basis for more complex model development. The long term stability of populations of 
single cell plants and animals (phytoplankton and zooplankton) is primarily based on two 
classical challenges. The first is the cleverness of the animals whose feeding or 
reproductive behavior produces density dependent responses which induce stability. The 
second challenge is that the complexity of the physical environment may play a critical 
role in preventing instabilities from growing either through dispersion or mixing of 
populations over a range of habitats. 

Steel (1977) derived the most basic equations for the concentrations of nutrients 
(N), phytoplankton (P) and herbivorous zooplankton (H). They are: 
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dN 

— — = -aNP + (6 - c)PH + k(N 0 - N) 
dt 



— = aNP - bPH - kP; 
dt 



dH 

— — = cPH - dH; 
dt 

where a is the phytoplankton growth rate and d is the zooplankton mortality rate. The 
grazing by H on P has the simple form of bPH, and the growth of H is some fraction of 
food intake ( c<b ). These equations are based on a simple two layered ocean (Figure 1) 
and the concentrations of N, P and H are normalized on standard or average values so that 
of N, P and H are dimensionless. The units of the model constants are in t' 1 . It is also 
assumed that all production takes place in the upper mixed layer since the lower layer is 
assumed to have a light intensity too low for production to occur and diurnal migrations 
of zooplankton are disregarded in these equations. This set of equations is more 
thoroughly described in Chapter HI and the model results are discussed in Chapter IV. 




Phytoplankton (P) 



Zooplankton (Z) 
diurnal migration 



P = 0 



Figure. 1 - A simple representation of a two-layer ecosystem. Reproduced from Steel 
( 1977 ) 
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Franks et. al (1986) developed a nutrient-phytoplankton- zooplankton (NPZ) 
model based on a modified version of the function experimentally determined by 
Mayzaud and Poulet in 1978. They compared their results to models based on the 
traditional Ivlev grazing expressions in which there is no asymptotic limit to the grazing 
rate. Their model showed highly damped oscillations in response to high grazing 
pressure as it approached steady state. 

In their conclusions they noted that in other experimental ecosystems such as 
CEPEX, MERL and Loch Ewe bag experiments that the biomass as a whole tended 
toward some equilibrium value. Species’ successions were seen in these experiments but 
the system did generally follow a non-oscillatory trend and approached some steady state. 
Franks believes that the Mayzaud and Poulet formulation allows a more predictable and 
robust model behavior and he suggested that it allows for a more realistic simulation of 
gross phytoplankton dynamics. Their results are displayed in Figure (2). 



5 




Time(doys) 



Figure 2. Time-dependent behavior of two types of herbivore grazing: (a) Ivlev and (b) 
Mayzaud and Poulet, from Franks et al. (1986). 



6 



Frost (1987) described seasonal variation of phytoplankton stock and production 
at ocean Station P (50° N 145° W). His model was designed to reproduce seasonal 
variations and integrated interactions of specific physical, chemical, and biological 
components of the subarctic pelagic ecosystem, and to simulate the observed level of 
phytoplankton production. 

In his model horizontal advection and diffusion were ignored. In addition, the 
phytoplankton were restricted to the mixed layer which was an oversimplification, but 
eliminated the need for detailed parameterization of mixing and phytoplankton growth 
and grazing mortality with in the pycnocline. He also assumed that, in general, 
phytoplankton are uniformly distributed in the mixed layer at Station “P” and there is no 
subsurface maximum; that is, the mixed layer is homogeneous with respect to 
phytoplankton. 

Fasham et al. (1990) presented a model of the annual cycles of plankton dynamics 
and nitrogen cycling in the oceanic mixed layer. This model was customized for Station 
“S” near Bermuda and utilized seven compartments or modules (Phytoplankton, 
Zooplankton, Bacteria, Nitrate, Ammonium, Dissolved organic nitrogen and Detritus). 
They found that one of the controlling parameters was the detritus sinking rates, for 
which they used values of 1 and 10 m/day. Due to the comparatively slow transition of 
detritus material into nutrients, it was found that the detritus only had a significant 
influence when sinking rates were on the order of 1 m/day. 

Fasham’s overall objective was to model the seasonal cycles of plankton and 
nutrients in the global ocean, in order to leam more about he role of oceanic organisms in 
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regulating the atmospheric CO 2 content. In his model, nitrogen is considered the limiting 
nutrient of primary production. The main physical forcing parameters in the model are 
the annual cycles of mixed layer depth and solar radiation. 

Fasham’s model is lacking in the physical aspects. His model does not include 
such processes as changes in turbulent mixing, effects of internal waves, seasonal 
changes in the diffusive mixing across the thermocline and within the mixed layer, and 
since it is only a 1-D model it does not include any horizontal advection. 

Aksnes and Lie (1990) presented a vertically resolved model for the land-locked 
fjord Lindaspollen, in western Norway. They found that primary production above the 
pycnocline is nutrient limited, while below the pycnocline the production is both nutrient 
and light limited. Hence, the new production during summer is controlled mainly by 
three factors: (1) the light availability in the deep nutrient-rich water, which depends on 
the intensity of the solar radiation and the turbidity of the water, (2) the turbulent 
diffusion of nutrients within the watermasses above the pycnocline, 
and (3) the nutrients supplied by freshwater runoff to the layer above the pycnocline. 

Aksnes and Egge (1991) discuss the difficulties in obtaining accurate nutrient 
uptake parameters and state that currently oceanographers are only trying to evaluate 
orders of magnitude vice actual rates. Nutrient uptake refers to processes necessary for 
moving a nutrient (denoted as an ion) from the outside to the inside of a cell (or perhaps 
to a site where further processing of the ion takes place.) The rate of nutrient uptake 
depends on the time period needed for transfer of nutrient ions from the outside to the 
inside of the cell, but also on the time elapsing between encounters of nutrient ions and 
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uptake sites. Increased temperature generally means higher biochemical rates. Thus it is 
reasonable that increased temperature will lead to shorter handling times because 
handling involves biochemical processes. Handling time is defined as the time interval 
where uptake of an ion at an uptake site is prevented because of handling of another ion. 

In a WHOI Technical Report, Davis and Steel (1994) reported on the findings of a 
five day workshop in June 1993 on the status of upper layer biological/physical modeling. 
A summary of critical issues was presented, and the following conclusions were drawn: 

• Horizontal variability is not well known 

• Lagrangian nature of flows is not well known 

• Future acoustic observational techniques will help resolve turbulence fields. 

During this workshop the members used three physical models, a simple 

conventional mixed layer model, PWP model Price et al. (1986) and the MY Mellor and 
Yamada (1982) level 2.5 model. These models were coupled with various versions of 
NPZ and NPZD biological models and some results of the various coupled models were 
presented. The report went on to discuss the effectiveness of the various mixed layer 
models. 

Stramska and Dickey (1994) present a simple short time scale NPZ model. Like 
other models their model cannot simulate effects related to horizontal variability, 
processes such as photoadaptation and photoinhibition, or biological diversity of 
phytoplankton populations. They also noted that because of the numerous physical and 
biological factors controlling both growth and losses of phytoplankton, models describing 
annual phytoplankton cycles must be rather exhaustive and contain considerable numbers 
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of equations and parameters. As a result, it is often difficult to evaluate to what degree 
the modeling results are affected by a particular assumption. In the case of extremely 
deep mixing events due to strong winds and surface heat loss, significant removal of 
phytoplankton stock from the euphotic zone occurs. 

Doney et al. (1996) designed a coupled 1-D biological/physical mixed layer model 
utilizing detailed boundary layer physics and relatively simple biological interactions. 
Doney noted that Fasham et al. (1990) used biological models that are too complex based 
on our limited knowledge of the coupling between ocean planetary boundary layer (PBL) 
physics and biological vertical structure and fluxes. Phytoplankton growth in Doney’s 
model is controlled by nutrient and light limitations, and phytoplankton losses occur via 
natural mortality, zooplankton grazing and phytoplankton aggregation. The latter two 
processes lead to the production of sinking detrital particles. Doney’s model incorporates 
four compartments for phytoplankton, zooplankton, nutrients and detritus. One 
difference that Doney applied was the variations in phytoplankton chlorophyll to nitrogen 
ratios (chl:N). 

A one dimensional coupled biological and physical model for the Black Sea was 
proposed and tested by Oguz et al. (1996). For the physical model they used the Mellor- 
Yamada level 2.5 closure scheme. The biological part of the model involved interactions 
between inorganic nitrogen (nitrates and ammonium), phytoplankton, herbivorous 
zooplankton and detritus. In this model nitrogen is considered to be the limiting factor in 
primary production. The phytoplankton production process is parameterized in terms of 
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Liebig's law which assumes either lack of nutrients or light controls the biological 
productivity. 

The biological model requires a choice of 16 externally imposed parameters. 

Most of these parameters are poorly known from the observations. Determination of the 
values of the rate constants is one of the most difficult problems which ongoing research 
and observations will continue to verify and improve in the future. 

Oguz found that the initiation of the spring bloom depends strongly on the local 
mixing conditions and follows the weakening of the convective overturning mechanism. 
As soon as the surface layer of the water column gains a slight stability, the bloom 
commences before the formation of the seasonal thermocline. This model is capable of 
reproducing basic features of the plankton-nutrient fields in the Black Sea and helps us 
understand and interpret the available observations. To make this model more realistic it 
would need the introduction of the high-frequency variability and thus incorporation of 
strong wind-induced mixing, and parameterization of intermittent lateral nutrient input 
into the surface layer. 

Smythe et al. (1997) quantified the decay rate of subsurface turbulence using data 
obtained during several squalls that occurred in the western equatorial Pacific during the 
westerly windburst of December 1992. The ocean was generally well-mixed down to 
about 70 m (the top of the main thermocline) due to persistent strong winds associated 
with the windburst. They found that a typical squall event coincided with the appearance 
of a layer of relatively cool, fresh, highly turbulent water in the upper few meters. The 
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layer remained distinct from the underlying water for several hours after the squall had 
passed, during which time it spread to a depth of 10-30 m. 

They took successive profiles of viscous dissipation rates (e) following each of 
the squalls and found that within the squall layer, dissipation rates increased by an order 
of magnitude coincident with the increase in surface forcing. Below the squall layer, an 
equally dramatic decrease in € is evident. 

The observed suppression of subsurface turbulence provided some clues about the 
manner by which the turbulence was originally generated. They hypothesized that the 
mechanism for local TKE maintenance was a downward flux of surface-generated 
turbulence. Prior to squall passage, weak stratification allowed surface-generated 
turbulence to be readily distributed over the upper few tens of meters. This suggests that 
at least a significant fraction of ambient turbulence was driven by wind and buoyancy 
forcing. However, with the onset of heavy precipitation, near-surface stability 
characteristics change dramatically. The base of the squall layer is characterized by 
intense static stability due to salinity stratification. This stable stratification acted to 
inhibit downward transport of turbulence, thereby concentrating surface-generated 
turbulence within the squall layer. With the flux of turbulence from the surface thus 
suppressed, turbulence levels in the region below the squall layer decayed. 
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III. MODEL DESCRIPTIONS 



The two coupled physical/biological models will be described in this section. 

Both of the models were written in MATLAB code. The physical half of the model 
remained the same while varying the complexity of the biological component. 

A. THE PHYSICAL MODEL 

The physical half of the model is solved using the Runge-Kutta method with a 
variable time step to solve the mixed layer equations. The numerical scheme maximizes 
the time step to solve the set of equations without exceeding specified error tolerances. In 
mixed layer physics when mixed layer shallowing occurs, 3h/dt is not a continuous 
function. Therefore, the Runge-Kutta solution can only be used during periods of mixed 
layer deepening. Due to the short time steps and the inclusion of TKE unsteadiness, this 
model is ideally suited for predicting transient events and spin-up and spin-down of 
turbulence. The basic physical differential equations used are: 

Mixed Layer Depth 




( 1 ) 



Total TKE 




( 2 ) 



Vertical TKE 




( 3 ) 
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Mixed Layer Temperature 



%T) _ Qo A f^e 
dt pc p h h 



( 4 ) 



Easterly Velocity 



d{hu) - t x 
— — = fhv + — 
a p 



(5) 



Northerly Velocity 



- T„ 



6 '(hv) - -y 

— — = -jhu + — 
a p 



( 6 ) 



Equation (1) defines the mixed layer depth, which can change due to the 
entrainment velocity, w e and any upwelling or down welling ( W z= - h ). Equations (2) and 
(3) calculate total turbulent kinetic energy (TKE) times h and vertical TKE times h 
respectively. Equation (4) calculates mixed layer temperature, which is governed by the 
surface heat flux and entrainment heat flux. Equations (5) and (6) are used to fmd the 
easterly velocity (u) times h and the northerly velocity (v) times h. 

For a more complete description of the physical processes refer to Stone (1997) 
and Garwood (1977). 



B. THE BIOLOGICAL MODEL 



The first version of the biological model is based on the equations developed by 
Steel (1977) based on the well-known Lotka-Volterra predator-prey relations. The 
equations have been modified to include terms for nutrient and phytoplankton 
entrainment. In this model, terms for Detritus (waste matter and decaying organisms 
which in turn produce nutrients) have been neglected. The equations are: 
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d 

— < N >= -a < P >< N > +{b - c) < P >< H > - 

dt 



d A Pw 

— <P>=a<PxN>-b<PxH> - — - — 

a h 



d 

— < H >= c < P x H > -d < H> 

a 



ANw e 

h 



(7) 

( 8 ) 
(9) 



where N is nutrients, P is phytoplankton, H is single-cell herbivores (or zooplankton), w e 
is entrainment and h is the mixed layer depth. The < > denotes the time average mean 
value. A conceptual view of this three component biological model is shown in Figure 3. 
In this simple model the phytoplankton consume the nutrients. Some of the 
phytoplankton are consumed by the zooplankton and others die and are returned to the 
system as nutrients. The zooplankton produce waste matter and eventually die, both by 
products are then returned to the system as nutrients. 




Figure 3. Schematic diagram of the three component biological model showing the flow 
pathways. 

Four constants are needed; the phytoplankton growth rate, (a), the growth by 
grazing term, (b), the recycled waste from grazing term, (c), and the zooplankton death 
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term (d). The values chosen for these were a = 2.0/day, b = 1.5/day, c = 1.2/day, and d = 
0.2/day. 

When determining how to use the nutrients, all of the elements that make up 
nutrients, such as nitrates, nitrites and phosphates, are grouped together in a general term 
as nutrients. This is an area where large generalizations are made for the biological 
model. I assumed a normalized value of 1 for concentration of nutrients below the mixed 
layer and a value less than 1 for concentration of nutrients within the mixed layer. To 
initialize the model, a value of 0.1 was chosen for the concentration of nutrients in the 
mixed layer. 

For the phytoplankton, the different species are combined into a single term. 
Since phytoplankton are very dependent on light they typically only reside in the euphotic 
zone (the area of greater than 1% light intensity of the surface light level). The mixed 
layer depth (MLD) and the euphotic zone are assumed to be the same. Therefore, the 
phytoplankton can only reside in the mixed layer (normalized concentration of 1) and are 
zero below the mixed layer. A value of 0.01 was chosen to initialize the model. 

Similar simplifications are made for modeling the zooplankton. Herbivorous 
zooplankton have the ability to propel themselves through the water and often have a 
diurnal cycle depending upon the species. Zooplankton are not dependent upon light and 
are not confined to the mixed layer. However, for the model, it is assumed that the 
zooplankton are not species specific and are going to remain in the vicinity of their food 
source (the phytoplankton) and are not going to have diurnal migrations. The 
concentration for zooplankton used to initialize the model is 0.01. 
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The second, more advanced biological model is based on the Davis, Olson and 
Flierl NPZD model outlined in Davis and Steel (1994) which is derived from the Flierl 
and Davis (1993) NPZ model. The set of coefficients was also taken from Davis and 
Steel (1994) and is outlined in Table 1. Most of the values are accepted for 
diatom/copepod dominated food webs, based on historical observations and experimental 
data. The coefficients with the lowest degree of accuracy are the zooplankton death rate 
(zd), the detritus remineralization rate (rd), and the phytoplankton death rate (pd). All 
other assumptions made for the previous biological model remained the same for this 
model. The schematic diagram of the model showing the pathways and directions of 
flows between the four compartments is displayed in Figure 2. 

In the second version of the model a solar radiation component was added and 
was modeled based on Frost (1987) and Oguz et al. (1996). This exponential relationship 
is: 

I(z) = I 0 exp(-(kw*(P* kc))z (10) 

where I(z) is the solar radiation at depth z. The model has the capacity to toggle I(z) on 
and off to simulate a diurnal variation. Io is the incident solar radiation at the surface 
corrected to only account for photosynthetically active radiation (PAR). There is some 
debate over the actual PAR value. Frost (1987) and Oguz et al. (1996) concluded the 
value to be 50%. Later Frost (1993) used a value of 70%, Asknes & Lie (1990) 
determined the value to be 25% while Skogen et al. (1995) used a PAR value of 40%. A 
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PAR of 50% was chosen for this study. The extinction coefficient for attenuation due to 
sea water is kw and kc*P represents the attenuation due to phytoplankton in the water 
column. Again there has been some debate over these coefficients. For this research, 
values of kw and kc were taken from Oguz et al. (1996) and are kw = 0.08 and kc = 0.07. 
A more thorough explanation of light attenuation characteristics in sea water may be 
found in White (1991) and Cullen (1990). 
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Figure 4. Schematic diagram of the four component biological model showing the flow 
pathways. Reprinted from Doney et al. (1996). 
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The advanced set of equations are as follows: 



Nutrients: 

dN 

= -pg* P* N * PAR + ( ae - ge ) *gr*P*Z + ae*zd*Z + rd*D- nent (11) 

dt 

Phytoplankton: 



— = pg* P* N * I - gr* P* Z - pd* P - pent 
dt 



( 12 ) 



Zooplankton: 



dZ 

= ge*gr*P*Z—zd*Z 



( 13 ) 



Detritus: 






( 14 ) 



where: 



nutrient uptake 
grazing 

phytoplankton death 
zooplankton death 
detritus remineralization 



= pg* P* N* PAR 
= gr*P*Z 
- pd* P 
= zd*Z 
= rd*D 



Table 1 defines the above terms, including units and initial conditions where 
applicable. 
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Biological Model Components and units 



DESCRIPTION 


SYMBOL 


INITIAL VALUE 


UNITS 


Nutrient Concentration (Summer) 


N 


0.1 


dimensionless | 


Nutrient Concentration (Winter) 


N 


1.0 


dimensionless 


Phytoplankton Concentration 


P 


0.01 


dimensionless 


Zooplankton Concentration 


Z 


0.01 


dimensionless 


Detritus Concentration 


D 


0.01 


dimensionless 


Solar Radiation (Summer) 


PAR 


100 


w/m A 2 


Solar Radiation (Winter) 


PAR 


75 


w/m A 2 


Phytoplankton Growth Rate 


Pg 


0.3 


/day 


Zooplankton Grazing Rate 


gr 


0.70 


/day 


Zooplankton Death Rate 


zd 


0.05 


/day 


Detritus Remineralization Rate 


rd 


0.2 


/day 


Phytoplankton Death Rate 


P d 


0.05 


/day 


Growth Efficiency for Z 


ge 


0.35 


dimensionless 


Assimilation efficiency for Z 


ae 


0.70 


dimensionless 



Table 1: Initial conditions for advanced biological model. 



In this model the phytoplankton growth rate is linearly dependent upon both light 
and nutrients. The zooplankton grazing rate is linearly dependent on the phytoplankton 
biomass. The detritus module is a simple “microbial loop” where the detritus 
remineralization rate (rd) is the rate at which bacteria convert detritus-nitrogen to 
dissolved nitrogen. 
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IV. RESULTS 



This study begins by using a simple NPZ model based on Steel (1977). The 
purpose of the first run was to compare to Franks et al. (1986) to the NPZ model used in 
this study. In this run an average mixed layer temperature was used vice a temperature 
“jump” at the bottom of the mixed layer. The zooplankton death rate (0.2/day) and the 
growth by grazing term (1.5/day) were taken from Franks et al. (1986). The other two 
initial condition variables were taken from an average of values given in several 
references. Figure 2 displays the results from the Franks model and Figure 5 displays the 
results from the NPZ model used here. 

The results of the two models are similar in that they both start with an initial 
oscillation of N, P and Z and then both approach a steady state. The Franks model attains 
a steady state in approximately 12 days as compared to approximately 20 days to attain a 
nearly steady state condition for the NPZ model used here. The differences in time are a 
result of the more complex biological interactions included in the Franks model. The key 
result here is the fact that both models do attain a quasi-steady state. 
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Mixed layer Biology - Comparison to Franks, et al. (1986) 




Figure 5. Mixed layer biology initialized in calm to moderate wind conditions. Q = 
0.001 cal/cm /s and t = 1.0 dynes/cm utilizing an average mixed layer temperature. For 
comparison to Franks et al. (1986). 
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In the next case the NPZ model was initialized in calmer conditions (wind stress 
(t) = 0.5 dynes/cm 2 and net heat flux (Q) = 0.001 cal/cm 2 /s). A highly stratified mixed 
layer was assumed with sea surface temperature of 20.0 °C and a -5.0 °C jump at the 
bottom of the mixed layer. The initial mixed layer depth was 1 0m. 

Figure 6a displays the evolution of the mixed layer biology from this run. Of 
interest is the fact that the system takes approximately 40 days to reach a nearly steady 
state condition. In addition, the concentration values are much lower than with a stronger 
wind stress. Figure 6b displays the TKE, mixed layer depth and horizontal currents for 
this case. In this case the mixed layer only reached a depth of 1 7 m in this 40 day run. 

The long time to attain a steady state, the low biological concentrations and the 
shallow mixed layer are a result of the low surface wind stress value. With a low wind 
stress the turbulent mixing is low and the entrainment of nutrients into the mixed layer is 
slow. As will be shown in the following figures, with increasing wind stress, the 
entrainment rapidly increases and the system will tend towards some steady state much 
more quickly. 
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Mixed layer Biology - Autumn - No storm - Low taux 




Figure 6a. Mixed layer biology initialized in calm conditions. Q - 0.001 cal/cm 2 /s and 
t = 0.5 dynes/cm 2 . 
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Figure 6b. TKE, mixed layer depth and horizontal currents initialized in calm conditions. 
Q = 0.001 cal/cm 2 /s and x = 0.5 dynes/cm 2 . 
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The remaining results are all based on a series of four sets of conditions imposed 
on each of the three model versions. The first being the same NPZ model without terms 
for light or detritus. The second version utilizes the NPZ equations with terms for light 
included and the final version utilizes a set of more advanced NPZD equations with a 
light budget and terms for detritus. This provides for a progression of complexity in the 
models studied here in order to draw some conclusion about the effects of light and 
detritus in the short term. 

The four sets of conditions are autumn with no storm induced during the run, 
autumn with a storm induce half way through the run, winter with no storm induced and 
winter with a storm induced. In each of the winter-time storm cases the mixed layer was 
found to deepen to very rapidly and to non-realistic depths. For that reason the winter 
storm cases are not presented here. Figure 7 graphically represents these four sets of 
conditions. 

In all cases the initial mixed layer depth is 10m and the system has a net 
downward heat flux of (Q) = 0.001 cal/cm /s. For the autumn conditions an initial sea 
surface temperature of 20 °C with a -5.0 °C temperature jump at the bottom of the mixed 
layer is utilized. The nutrient concentration within the mixed layer is initially 0.1 and 
maximum of 1 .0 below the mixed layer. 

In the winter cases a much weaker stratification exists with the sea surface 
temperature of 1 5. 1 °C and a temperature jump of 0. 1 °C at the bottom of the mixed layer. 
In this winter case the initial nutrient level is 1 .0 within the mixed layer. 
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In the cases where a storm is induced, the storm forcing is immediately at a 

'y 

maximum value of wind stress (r) = 2.0 dynes/cm and net upward heat flux ( Q ) = -0.006 
cal/cm 2 /s. The storm condition is utilized for two days and then the wind stress value 

'S 

returns to (r) = 1.0 dynes/cm ; however the net heat flux is kept in the upward direction to 
prevent the model from trying to shallow the mixed layer depth, which would cause 
inaccurate results. 
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Figure 7. Schematic diagram of the autumn and winter initial conditions and the storm 
conditions induced for these experiments. 
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The next four cases utilize the NPZ model with no terms for light or detritus. The 
initial mixed layer depth is 10m in all cases. Figure 8a shows the biological response 
with simulated autumn conditions of wind stress (r) = 1.0 dynes/cm 2 and net downward 
heat flux (0 = 0.001 cal/cm 2 /s. No storm was induced in this run. 

This case represents a standard baseline case for the next seven runs which all are 
some variation of this run. Note that the nutrient peak occurs at approximately day 4 and 
the phytoplankton peak at day 10, followed by the zooplankton peak at day 1 7. 

Maximum concentrations are approximately 0.43, 0.62 and 0.33 respectively. The system 
approaches a steady state at day 25. 

Figure 8b displays the TKE, mixed layer depth and horizontal currents. Note that 
the mixed layer deepened to 70 m in this case. 
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Mixed layer Biology - Autumn - No storm 




Figure 8a. Mixed layer biology initialized in autumn calm conditions (Q = 0.001 
cal/cm 2 /s and x = 1 .0 dynes/cm 2 ) and a strongly stratified water column. 
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Figure 8b. TKE, mixed layer depth and horizontal currents initialized in autumn calm 
conditions (Q = 0.001 cal/cm /s and x = 1 .0 dynes/cm ) and a strongly stratified water 
column. 
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The next run was initialized using the same conditions. In this run, however, a 
storm was induced on day 20 and stopped on day 22. The storm consisted on an 
increased wind stress (t) = 2.0 dynes/cm 2 and a strong upward heat flux of (Q) = -0.006 
cal/cm 2 /s. After the storm was stopped on day 22 the wind stress was returned to pre- 
storm conditions and the heat flux was kept in the upward direction but with a magnitude 
of (Q) = -0.001 cal/cm 2 /s. 

Figure 9a displays the biological conditions through this 40 day run. When the 
storm is induced on day 20 there is sharp increase in the nutrient level as a result of the 
extra nutrients being entrained from below the mixed layer. Likewise, the phytoplankton 
and zooplankton peak four and nine days later respectively. Of interest in this Figure is 
the fact that the phytoplankton peak is lower than the zooplankton peak in the post storm 
conditions. This is due to increased level of zooplankton when the storm is initialized 
which consume the phytoplankton more rapidly than prior to the storm. 

Figure 9b displays a dramatic increase in TKE when the storm is initiated at day 
20. The mixed layer deepens slowly prior to the storm and then deepens rapidly after the 
storm starts. Since there is no resistance to the deepening once it has started, the mixed 
layer continues to deepen rapidly to a depth of 400 m. There is also a perturbation in the 
horizontal currents once the storm is initiated. 
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Mixed layer Biology - Autumn - Storm 




days 



Figure 9a. Mixed layer biology initialized in calm autumn conditions(Q = 0.001 
cal/cm /s and x = 1 .0 dynes/cm ) with storm induced on day 20 (Q = -0.006 cal/cm /s and 
x = 2.0 dynes/cm 2 ) and a strongly stratified water column. 
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Figure 9b. TKE, mixed layer depth and horizontal currents initialized in calm autumn 
conditions(Q = 0.001 cal/cm 2 /s and x = 1.0 dynes/cm 2 ) with storm induced on day 20 (Q 
= -0.006 cal/cm 2 /s and x = 2.0 dynes/cm 2 ) and a strongly stratified water column. 
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The next case is identical to the initial conditions that produced the results in 
Figures 8a and 8b except that these cases are initialized using a weakly stratified water 
column to represent winter-time conditions. This was done in the model by using a 
temperature jump of -0.1 °C at the bottom of the mixed layer. The level of nutrients was 
initialized with a maximum value of 1 .0. All other initial conditions were the same. 

Comparing Figures 8a and 10a one can easily see that the system approaches a 
steady state much more quickly. The phytoplankton peak in Figure 8a is at day 10 where 
as in Figure 1 0a the phytoplankton peak occurs at day six. Likewise the zooplankton 
peak occurs at day 1 7 in Figure 8a vice day 1 1 in Figure 1 0a. 

Figure 10b shows the mixed layer deepening from the initial 10 m to 300 m in 40 
days. Obviously, with this deep of a mixed layer light limiting the phytoplankton growth 
plays a very important role in attempting to accurately predict the biological reaction to 
these conditions. The attenuation of light would prohibit the growth of phytoplankton at 
only a few tens of meters below the surface. For this reason this case it believed to be 
grossly inaccurate and is only shown here to justify the necessity of including a light 
budget into the model. A winter time case with a storm induced will not be discussed 
here for the reason mentioned above. 
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Mixed layer Biology - Winter - No storm 




Figure 10a. Mixed layer biology initialized in autumn calm conditions (Q = 0.001 
cal/cm2/s and t = 1 .0 dynes/cm 2 ) and a weakly stratified water column. 
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Figure 10b. TKE, mixed layer depth and horizontal currents initialized in calm autumn 
conditions(Q = 0.001 cal/cm 2 /s and t = 1.0 dynes/cm 2 ) and a weakly stratified water 
column. 
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In the next cases a light budget is added to the previously described NPZ model. 
The initial conditions are wind stress t = 1 .0 dynes/cm , a net downward heat flux of Q = 
0.001 cal/cm 2 /s, a temperature jump of 5 °C at the bottom of the mixed layer and initial 
nutrient concentration of 0. 1 . Figure 11a displays the evolution of the mixed layer 
biology. The phytoplankton peak occurs at day 5 with a concentration of approximately 
0.53. Compared to Figure 8a, the phytoplankton peak occurs 5 days earlier at a 20% 
higher concentration. 

The nutrient level in this case drops to nearly zero very rapidly after the start of 
the run as a result of rapidly increasing phytoplankton which consume the nutrients. The 
zooplankton peak occurs with the same time interval from the phytoplankton peak as in 
the case with no light budget. This is because the zooplankton have no response to light 
conditions in this model. 

As would be expected the TKE, mixed layer depth, and horizontal currents in 
Figure 1 lb are the same as in the model without the light budget (Figure 8b). 
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Mixed layer Biology - Autumn - No storm - light budget included 




days 



Figure 1 la. Mixed layer biology initialized in autumn calm conditions (Q = 0.001 
cal/cm 2 /s and x = 1 .0 dynes/cm 2 ) and a strongly stratified water column with a light 
budget included. 
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Figure 1 lb. TKE, mixed layer depth and horizontal currents initialized in autumn calm 
conditions (Q = 0.001 cal/cm 2 /s and x = 0.5 dynes/cm 2 ) and a strongly stratified water 
column with light budget included. 
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The next case was initialized with the same initial conditions as in the previous 
run. On day 20 a storm was initiated with a wind stress (x) = 2.0 dynes/cm 2 and net 
upward heat flux of (Q) = -0.006 cal/cm 2 /s. On day 22 the wind stress was returned to 
pre-storm conditions and the heat flux reduced to (Q) = -0.001 cal/cm 2 /s in the upward 
direction. 

Figure 12a shows the nutrient level growing rapidly when the storm is initiated 
because of the entrainment of nutrients from below the mixed layer as a result of 
increased turbulence. Figure 12b show the mixed layer deepening rapidly after day 20 to 
a maximum depth of 400 m in 40 days. As a result of this deep mixed layer the light 
cannot provide the PAR to maintain the growth of phytoplankton and subsequently the 
zooplankton, which both approach zero. 

This case demonstrates the definite necessity of a light budget to accurately model 
biological components. 
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Concentrations 



Mixed layer Biology - Autumn - Storm - light budget included 




Figure 12a. Mixed layer biology initialized in calm autumn conditions(Q = 0.001 

7 9 7 

cal/cm /s and x = 1 .0 dynes/cm ) with storm induced on day 20 (Q = -0.006 cal/cm /s and 

7 

x = 2.0 dynes/cm ) and a strongly stratified water column including a light budget. 
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Figure 12b. TKE, mixed layer depth and horizontal currents initialized in calm autumn 
conditions(Q = 0.001 cal/cm 2 /s and t = 1.0 dynes/cm 2 ) with storm induced on day 20 (Q 
= -0.006 cal/cm /s and t = 2.0 dynes/cm ) and a strongly stratified water column 
including a light budget. 
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The next case represents winter conditions with a light budget included in the 
NPZ model. The model was initialized using downward heat flux (Q) = 0.001 cal/cm 2 /s 
and wind stress (x) = 1 .0 dynes/cnT. The initial nutrient concentration was 1 .0. The water 
column was weakly stratified with temperature jump of 0.1 °C at the base of the mixed 
layer. 

Although the biological results in Figure 13a appear to be unrealistic, the response 
of the model is accurate. Figure 13b shows the mixed layer deepening rapidly which 
exceeds 1 00 m in less than three days. As a result of this rapidly deepening mixed layer 
the growth of phytoplankton is decreased and an increase in the level of nutrients occurs 
quickly after an initial consumption from the only growth of phytoplankton in this run. 
The zooplankton never get the chance to grow since the phytoplankton do not achieve a 
very high level. Also of note in Figure 13b is the rapid damping of the horizontal 
currents, which is a result of the rapidly deepening mixed layer. 

This case was rerun with a storm condition induced at day 20; however, the results 
are not shown since mixed layer deepens rapidly as in the previous run and when the 
storm was initialized it deepened to unrealistic depths. 
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Mixed layer Biology - Winter - No storm 




Figure 13a. Mixed layer biology initialized in winter conditions (Q = 0.001 cal/cm 2 /s 
and x = 1 .0 dynes/cm ) and a weakly stratified water column with a light budget included 
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Figure 13b. TKE, mixed layer depth and horizontal currents initialized in winter 
conditions (Q = 0.001 cal/cm 2 /s and x = 1 .0 dynes/cm 2 ) and a weakly stratified water 
column with light budget included. 
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In the next series of runs the NPZD model was used. This model included terms 



for detritus and light. The first of these runs is the basic autumn previously mentioned 
with (Q) = 0.001 cal/cm /s and (x) = 1.0 dynes/cm . The initial nutrient level is 0.1 and 
the water column is strongly stratified at the base mixed layer with a temperature jump of 
5.0 °C. The NPZD runs are conducted for 12 days vice 40 days for the NPZ model since 
the NPZD model approaches a steady state much more quickly than the NPZ model. 

Figure 14a displays the evolution of the nutrients, phytoplankton, zooplankton and 
detritus for this NPZD model run. The nutrient level peaks followed by a phytoplankton 
peak. As a result of this phytoplankton peak, the nutrient level shows an associated dip 
from the increased phytoplankton grazing. The phytoplankton concentration then dips at 
day five as a result of the growing zooplankton. The detritus concentration increases to a 
value of just greater than 0.05 in five days. The entire system approaches a steady state 
in approximately seven days. 

Figure 14b displays the TKE, mixed layer depth and horizontal currents for this 
12 day run. The reader will note from the figure that the mixed layer deepens from 10 m 
to 32 m in 12 days, so light limiting is not a major factor in this case. 
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Mixed layer Biology - Autumn - No Storm 




days 



Figure 14a. NPZD mixed layer biology initialized in autumn calm conditions (Q = 0.001 
cal/cm 2 /s and x = 1.0 dynes/cm 2 ) and a strongly stratified water column with a light 
budget included. 
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Figure 14b. NPZD TKE, mixed layer depth and horizontal currents initialized in autumn 
calm conditions (Q = 0.001 cal/cm 2 /s and x = 1.0 dynes/cm 2 ) and a strongly stratified 
water column with light budget included. 
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The next case is initialized with the same initial conditions as in the last case. For 



•y 'y 

this run a storm condition (Q) = -0.006 cal/cm Is and (x) = 2.0 dynes/cm started on day 
five and stopped on day seven. The post-storm conditions were (Q) = -0.001 cal/cm 2 /s 
and (x) = 1 .0 dynes/cm 2 . 

Figure 15a displays the NPZD mixed layer biology for this storm case. The most 
striking feature in this plot is the rapid increase in nutrients when the storm is initialized. 
Also of interest is the rapid decay of phytoplankton, zooplankton and detritus after the 
storm is started. An investigation of Figure 15b reveals the reason for this response. The 
mixed layer, deepening slowly for the first five days, begins to rapidly deepen after the 
storm initialization. As was shown in the previous rims with the NPZ model with the 
light budget, when the mixed layer deepens beyond approximately 100 m the 
phytoplankton growth is essentially cut off. 

Figure 15b also reveals a dramatic jump in TKE and a strong damping of 
horizontal currents once the storm condition is started as a result of the rapidly deepening 
mixed layer. 
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Mixed layer Biology - Autumn - Storm - NPZD 




Figure 15a. NPZD mixed layer biology initialized in calm autumn conditions(Q = 0.001 
cal/cm 2 /s and x = 1.0 dynes/cm 2 ) with storm induced on day five (Q = -0.006 cal/cm 2 /s 
and x = 2.0 dynes/cm 2 ) and a strongly stratified water column including a light budget. 
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Figure 15b. NPZD TKE, mixed layer depth and horizontal currents initialized in calm 
autumn conditions(Q = 0.001 cal/cm 2 /s and x = 1.0 dynes/cm 2 ) with storm induced on day 
five (Q = -0.006 cal/cm 2 /s and t = 2.0 dynes/cm 2 ) and a strongly stratified water column 
including a light budget. 
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For the final case, the NPZD model was forced with winter conditions Q = 0.001 
cal/cm /s and t = 1.0 dynes/cm . The water column is weakly stratified with a 
temperature jump of 0.1 °C and the initial nutrient concentration level of 1.0. 

Figure 16a displays the NPZD mixed layer biology for this winter time case. This 
plot indicates a very rapid approach to a steady state condition. As in previous cases the 
cause for this rapid change is the quickly deepening mixed layer as is displayed in Figure 
1 6b. The phytoplankton concentration instantaneously grows to a maximum value and 
then decays very rapidly on the order of one day. Similarly the nutrient concentration 
profile is the inverse of the phytoplankton profile. The zooplankton and detritus also 
peak rapidly with much lower concentration values compared to the phytoplankton. 

Figure 16b shows the TKE, mixed layer depth and horizontal currents for this run. 
Of interest is the rapid dampening of the horizontal currents as the mixed layer deepens 
rapidly. 

This case was also run with a storm condition induced; however, when the storm 
was induced the mixed layer again deepened very rapidly to an unrealistic depth. For this 
reason the results of that run are not described here. 
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Mixed layer Biology - Winter - No storm - NPZD 




Figure 16a. NPZD mixed layer biology initialized in winter conditions (Q = 0.001 
cal/cm /s and t = 1 .0 dynes/cm ) and a weakly stratified water column with a light budget 
included. 
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Figure 16b. NPZD TKE, mixed layer depth and horizontal currents initialized in winter 
conditions (Q = 0.001 cal/cm /s and x = 1.0 dynes/cm ) and a weakly stratified water 
column with light budget included. 
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V. CONCLUSIONS AND AREAS FOR FURTHER RESEARCH 



Many conclusions may be drawn from the results presented in the previous 
chapter. The most dramatic result which was seen throughout the model simulations was 
the fact that nutrient concentration levels increase rapidly as turbulence spins up because 
of the entrainment of the nutrients from below the mixed layer. 

It is also evident that stronger turbulence promotes a rapid approach to steady 
state conditions. This is result of the additional nutrients being introduced to the system 
and the subsequent rapid increase of phytoplankton and zooplankton. Surprisingly, the 
increasing nutrients drive the system to a near steady state more rapidly but do not 
increase the overall concentration levels of nutrients, phytoplankton or zooplankton. 

The light limiting the growth of phytoplankton is the major limiting factor below 
100 m in the experiments run here. This depth would depend greatly upon the clarity of 
the water mass being studied. It is evident from these results that future coupled/physical 
biological modeling must include a light budget. 

From the results presented here, it is obvious that turbulence induces a biological 
response. This response could be one cause of the initiation of a bloom. This finding is 
in concurrence with the findings of Oguz et al. ( 1 996) in which they found that the onset 
of a bloom in the Black Sea was, at least, partially dependent upon local mixing. 

As stated in Doney et al. (1996) our knowledge of planetary boundary layer 
physics and biological vertical structure and fluxes is limited. This fact dictates that the 
future of physical/biological modeling requires the improvement of the physical half of 
the model before we can attempt to make major improvements to the biological half. As 
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such, for advancements to take place in this field we must improve our knowledge of 
turbulence fields and our ability to accurately model them. As our knowledge of the 
complex interactions between biological organisms improves, so will our ability to build 
the groundwork for improved biological modeling. 

The physical/biological modeling field will be improved with the development of 
a high resolution in-situ data set, including detailed measurements of biological and 
physical parameters through a transient event, much like the physical data set collected by 
Smythe et al. (1997). This would be a great step in the advancement of model tuning and 
verification. 
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